Análise de Focos de Incêndio na Amazônia Legal

Sobre os dados

Focos de queimada

Os dados de focos foram obtidos do banco de dados do INPE: https://queimadas.dgi.inpe.br/queimadas/bdqueimadas#exportar-dados

Os dados foram filtrados pelos limites Amazônia legal, entre os anos de 2011 e 2020.

Outras variáveis

A distância à rodovia mais próxima foi calculada em função de um shapefile de rodovias disponibilizado pelo IBGE. O mapa de grau de urbanização também foi obtido do IBGE. Os dados de áreas desmatadas são elaborados pelo INPE através do PRODES.

Imports

Os pacotes sf e terra são essenciais para a manipulação de dados vetoriais e matriciais respectivamente. O pacote tidyverse contém funções uteis para a manipulação de tabelas. O pacote spatstat é utilizado para extrair estatisticas espaciais como a densidade de focos. O pacote caret é responsável pelo treinamento e aplicação dos modelos.

library(terra) 
library(tidyverse)
library(sf)
library(spatstat)
library(caret)
library(doParallel)
library(gbm)

Criação da camada focos/km²

A ideia geral aqui é criar uma imagem raster com a contagem da ocorrência de focos por km². Primeiro os dados de foco devem ser carregados. Eu aproveito para ler o shapefile da amazônia legal também.

focos <- st_read("focos_2011_2020.shp")

## Reading layer `focos_2011_2020' from data source 
##   `F:\Pablo\Documentos\Doutorado\side\focos_2011_2020.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1215884 features and 15 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 2828897 ymin: 8004499 xmax: 6111664 ymax: 10573720
## Projected CRS: SIRGAS 2000 / Brazil Polyconic

am_legal <- st_read("amazonia_legal.shp")

## Reading layer `amazonia_legal' from data source 
##   `F:\Pablo\Documentos\Doutorado\side\amazonia_legal.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 2794537 ymin: 8004219 xmax: 6112012 ymax: 10586380
## Projected CRS: SIRGAS 2000 / Brazil Polyconic

am_vect <- vect(am_legal)

Em seguida criar um raster vazio de resolução de 10km para preencher os dados:

focos_km2 <- rast(ext=ext(focos), crs=crs(focos), resolution=10000)

E enfim contar os focos por pixel da imagem:

focos_km2 <- rasterize(vect(focos), focos_km2, fun=length, field=1, background=0) %>%
  mask(am_vect)
colorPal <- colorRampPalette(rev(RColorBrewer::brewer.pal(11, "RdYlGn")))(256)
plot(main="Focos brutos por 10km²", focos_km2, col=colorPal)

Ou então fazer uma análise de densidade:

pts <- st_coordinates(focos) %>%
  unique()
pts_ppp <- ppp(pts[,1], pts[,2], window=as.owin(am_legal), check=FALSE)

A densidade é diferente da quantidade bruta de focos por pixel no sentido de que ela é suavizada em função de uma banda (distância). A função padrão de suavização é a gaussiana.

Um teste rapido usando diferentes bandas:

bandas <- c(10000, 25000, 50000, 100000)
dens <- lapply(bandas, function(x) rast(density(pts_ppp, sigma=x))*1e6)
dens <- rast(dens)
names(dens) <- c("10km", "25km", "50km", "100km")
plot(dens, col=colorPal)

Nesse caso uma banda de 50km pareceu apropriada.

pt_dens <- density(pts_ppp, sigma=50000, eps=10000)
pt_dens <- rast(pt_dens)*1e6
crs(pt_dens) <- crs(focos_km2)
pt_dens <- resample(pt_dens, focos_km2)
plot(main="Densidade de focos por km²", pt_dens, col=colorPal)

Lado a lado:

r_stack <- rast(list(focos_km2, pt_dens))
names(r_stack) <- c("Focos por 10km²", "Densidade de focos")
plot(r_stack, col=colorPal)

É importante lembrar que os valores calculados representam um histórico de 10 anos, portanto o valor anual seria a média, ou o valor apresentado dividido por 10.

Finalmente, salvar os resultados:

writeRaster(focos_km2, "focos_10km.tif", gdal=c("COMPRESS=DEFLATE", "PREDICTOR=2"), overwrite=T)
writeRaster(pt_dens, "densidade.tif", gdal=c("COMPRESS=DEFLATE", "PREDICTOR=2"), overwrite=T)

Variáveis independentes

Os mapas dessas variáveis ja foi feito usando o QGIS, então é só carregar eles.

variaveis <- list(
  grau_urb = rast("grau_urb.tif"),
  dist_rod = rast("dist_rod.tif"),
  dist_desm = rast("dist_desm.tif"),
  dist_ind = rast("dist_ind.tif"),
  dist_cons = rast("dist_cons.tif")
)

variaveis <- lapply(variaveis, function(x){
  resample(x, variaveis[[1]]) %>%
    mask(am_vect)
}) 

variaveis <- rast(variaveis)
plot(variaveis, col=viridis::viridis(256))

Amostragem

A princípio a amostragem vai ser feita por via de pontos aleatórios. Por enquanto vamos usar 5000 pontos:

set.seed(1)
amostras <- st_sample(am_legal, 5000) %>%
  st_sf('ID' = seq(length(.)), 'geometry' = .)
ggplot()+
  geom_sf(data=am_legal, col="gray", fill=NA)+
  geom_sf(data=amostras, col="red", pch= 3, alpha=0.3)+
  theme_bw()

Extração dos dados rasteriais e conversão para tabela:

amostras <- vect(amostras)
dados_y <- tibble(
  focos_10km2=terra::extract(focos_km2, amostras)[,2],
  densidade=terra::extract(pt_dens, amostras)[,2]
)
dados_x <- terra::extract(variaveis, amostras) %>%
  as_tibble() %>%
  select(-ID)
dados <- bind_cols(dados_y,dados_x)

Traduzindo em forma de gráfico:

dados %>%
  select(-focos_10km2) %>%
  mutate(across(-c(densidade, grau_urb), ~ .x/1000)) %>%
  pivot_longer(-densidade) %>%
  ggplot(aes(value, densidade))+
    geom_point(alpha=0.2)+
    xlab("Valor")+
    ylab("Densidade")+
    facet_wrap(vars(name), ncol = 3, scales = "free_x")

Alguns parâmetros pra tunagem dos modelos

fitControl <- trainControl(## 10-fold CV
                           method = "repeatedcv",
                           number = 10,
                           ## repeated ten times
                           repeats = 10)

rfGrid <- expand.grid(mtry = c(2,3,4,5))

gbmGrid <-  expand.grid(interaction.depth = c(1, 5, 9), 
                        n.trees = (1:30)*50, 
                        shrinkage = 0.1,
                        n.minobsinnode = 20)

Regressão com RandomForest

cl <- makePSOCKcluster(5)
registerDoParallel(cl)
modelo_rf <- dados %>% 
  drop_na(densidade) %>%
  train(densidade ~ grau_urb + dist_rod + dist_desm + dist_ind + dist_cons, data=.,
        trControl = fitControl,
        tuneGrid = rfGrid,)
stopCluster(cl)
write_rds(modelo_rf, "modelo_rf.rds")

Estatisticas do modelo:

modelo_rf

## Random Forest 
## 
## 4955 samples
##    5 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 10 times) 
## Summary of sample sizes: 4459, 4459, 4459, 4461, 4459, 4459, ... 
## Resampling results across tuning parameters:
## 
##   mtry  RMSE       Rsquared   MAE       
##   2     0.1270520  0.6769336  0.08797339
##   3     0.1255838  0.6830307  0.08581174
##   4     0.1248307  0.6863201  0.08456322
##   5     0.1243623  0.6884374  0.08367233
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 5.

plot(modelo_rf)

plot(varImp(modelo_rf))

Comparação com a realidade:

dados2 <- dados %>%
  mutate(preds = predict(modelo_rf, newdata=.), erro=(densidade-preds))
dados2 %>%
  select(-focos_10km2) %>%
  pivot_longer(-c(densidade, preds, erro)) %>%
  pivot_longer(-c(name, erro, value), names_to = "tipo", values_to = "densidade") %>%
  ggplot(aes(value, densidade, color=tipo))+
    geom_point(alpha=0.3)+
    facet_wrap(vars(name), ncol=2, scales="free_x")+
    scale_color_discrete("Legenda: ",labels=c("Previsto", "Observado"))+
    theme(legend.position = "bottom")

Distribuição de erros:

dados2 %>%
  ggplot(aes(x=erro))+
  geom_histogram(bins=50)+
  scale_x_continuous(breaks=seq(-1,1,0.1))

Enfim classificar a imagem:

pred_rast <- predict(variaveis, modelo_rf, na.rm=T, cores=2)
writeRaster(pred_rast, "rf_rast.tif", gdal=c("COMPRESS=DEFLATE", "PREDICTOR=2"), overwrite=TRUE)

plot(pred_rast, col=colorPal)

par(mfrow=c(1,2))
plot(main="Previsto", pred_rast, col=colorPal)
plot(main="Observado", pt_dens, col=colorPal)

Regressão com Stochastic Gradient Boosting (GBM)

cl <- makePSOCKcluster(5)
registerDoParallel(cl)
modelo_gbm <- dados %>% 
  drop_na(densidade) %>%
  train(densidade ~ grau_urb + dist_rod + dist_desm + dist_ind + dist_cons,
        data=.,
        method = "gbm",
        trControl = fitControl,
        tuneGrid = gbmGrid,
        verbose = FALSE)
stopCluster(cl)
write_rds(modelo_gbm, "modelo_gbm.rds")

Estatisticas do modelo:

plot(modelo_gbm)

plot(varImp(modelo_gbm))

Comparação com a realidade:

dados3 <- dados %>%
  mutate(preds = predict(modelo_gbm, newdata=.), erro=(densidade-preds))
dados3 %>%
  select(-focos_10km2) %>%
  pivot_longer(-c(densidade, preds, erro)) %>%
  pivot_longer(-c(name, erro, value), names_to = "tipo", values_to = "densidade") %>%
  ggplot(aes(value, densidade, color=tipo))+
    geom_point(alpha=0.3)+
    facet_wrap(vars(name), ncol=2, scales="free_x")+
    scale_color_discrete("Legenda: ",labels=c("Previsto", "Observado"))+
    theme(legend.position = "bottom")

Distribuição de erros:

dados3 %>%
  ggplot(aes(x=erro))+
  geom_histogram(bins=50)+
  scale_x_continuous(breaks=seq(-1,1,0.1))

Aplicação do modelo:

pred_rast <- predict(variaveis, modelo_gbm, na.rm=T, cores=2)
writeRaster(pred_rast, "gbm_rast.tif", gdal=c("COMPRESS=DEFLATE", "PREDICTOR=2"), overwrite=TRUE)

plot(pred_rast, col=colorPal)

par(mfrow=c(1,2))
plot(main="Previsto", pred_rast, col=colorPal)
plot(main="Observado", pt_dens, col=colorPal)

Extreme gradient boost

cl <- makePSOCKcluster(5)
registerDoParallel(cl)
modelo_xgb <- dados %>% 
  drop_na(densidade) %>%
  train(densidade ~ grau_urb + dist_rod + dist_desm + dist_ind + dist_cons,
        data=.,
        method = 'xgbTree',
        trControl = fitControl,
        verbose = FALSE)
stopCluster(cl)
write_rds(modelo_xgb, "modelo_xgb.rds")

plot(modelo_xgb)

dados4 <- dados %>%
  mutate(preds = predict(modelo_xgb, newdata=.), erro=(densidade-preds))

dados4 %>%
  select(-focos_10km2) %>%
  pivot_longer(-c(densidade, preds, erro)) %>%
  pivot_longer(-c(name, erro, value), names_to = "tipo", values_to = "densidade") %>%
  ggplot(aes(value, densidade, color=tipo))+
    geom_point(alpha=0.3)+
    facet_wrap(vars(name), ncol=2, scales="free_x")+
    scale_color_discrete("Legenda: ",labels=c("Previsto", "Observado"))+
    theme(legend.position = "bottom")

Distribuição de erros:

dados4 %>%
  ggplot(aes(x=erro))+
  geom_histogram(bins=50)+
  scale_x_continuous(breaks=seq(-1,1,0.1))

Aplicação do modelo:

pred_rast <- predict(variaveis, modelo_xgb, na.rm=T, cores=2)
writeRaster(pred_rast, "xgb_rast.tif", gdal=c("COMPRESS=DEFLATE", "PREDICTOR=2"), overwrite=TRUE)

plot(pred_rast, col=colorPal)

Comparação entre os erros dos modelos:

ggplot()+
  geom_freqpoly(aes(dados2$erro, color="RF"), bins=100)+
  geom_freqpoly(aes(dados3$erro, color="GBM"), bins=100)+
  geom_freqpoly(aes(dados4$erro, color="XGB"), bins=100)+
  theme_bw()+
  xlab("Erro")+
  ylab("Frequência")+
  theme(legend.title = element_blank())

rmse_modelos <- tibble(
  "RF" = min(modelo_rf$results$RMSE),
  "GBM" = min(modelo_gbm$results$RMSE),
  "XGB" = min(modelo_xgb$results$RMSE)
) %>%
  pivot_longer(everything())

rmse_modelos %>%
  arrange(value) %>%
  ggplot(aes(reorder(name, -value), value))+
    geom_col(fill="black", width=0.3)+
    ylab("RMSE")+
    xlab("Modelo")